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Abstract — We study Voronoi diagrams for distance functions 
that add together two convex functions, each taking as its argu- 
ment the difference between Cartesian coordinates of two planar 
points. When the functions do not grow too quickly, then the 
Voronoi diagram has linear complexity and can be constructed 
in near-linear randomized expected time. Additionally, the level 
sets of the distances from the sites form a family of pseudocircles 
in the plane, all cells in the Voronoi diagram are connected, 
and the set of bisectors separating any one cell in the diagram 
from each of the others forms an arrangement of pseudolines in 
the plane. We apply these results to the smoothed distance or 
biotope transform metric, a geometric analogue of the Jaccard 
distance whose Voronoi diagrams can be used to determine the 
dilation of a star network with a given hub. For sufficiently closely 
spaced points in the plane, the Voronoi diagram of smoothed 
distance has linear complexity and can be computed efficiently. 
We also experiment with a variant of Lloyd's algorithm, adapted 
to smoothed distance, to find uniformly spaced point samples with 
exponentially decreasing density around a given point. 

Index Terms — biotope transform metric; convex function; di- 
lation; Lloyd's algorithm; pseudocircle; pseudoline; randomized 
algorithms smoothed distance; Voronoi diagram. 

I. Introduction 

Any bivariate function / and finite set of point sites p. L in 
the plane give rise to a minimization diagram in which the 
cell for site pi consists of points q such that the value of the 
translated function f(q — pi) is less than or equal to the value 
of any of the other translates of /. A familiar example is given 
by Euclidean Voronoi diagrams, the minimization diagrams of 
translates of the convex functions f{x,y) — y/ x 2 + y 2 (or, 
equivalently, f(x,y) — x 2 + y 2 ) that measure the (squared) 
distance of (x, y) from the origin. 

Minimization diagrams have many applications, and typi- 
cally have quadratic complexity |9), but in some important 
cases their complexity is much smaller. In a Euclidean Voronoi 
diagram, each cell is a convex polygon, so the Voronoi 
diagram for n sites partitions the plane into n connected 
regions and has 0{n) vertices and edges. These combinatorial 
facts form the basis of efficient algorithms for constructing 
Euclidean Voronoi diagrams and using them in other geometric 
algorithms. It is natural, therefore, to ask: Which other mini- 



mization diagrams have connected cells and a linear number 
of number of features? 

A partial answer was provided by Chew and Drysdale (Tl: 
Voronoi diagrams for convex distance functions have con- 
nected cells and linear complexity. A convex function / is a 
convex distance function if, for any positive scalar a and any 
point p, f(ap) — af(p). Not every convex function has this 
property: it implies that the level sets {p \ f(p) — k} are all 
similar, and that the function is linear on rays from the origin, 
neither of which are true for all convex functions. Another 
answer is given by the abstract Voronoi diagrams p4| defined 
by a family of bisector curves that are required to intersect 
each other finitely many times and form simply connected 
cells. Bregman Voronoi diagrams pO) fall into this class: 
they have linear bisectors and convex-polygon cells. Abstract 
Voronoi diagrams may be constructed efficiently [ |T4| , fl5) , 
1 19] but it is unclear how to tell whether a given convex 
function has minimization diagrams that form abstract Voronoi 
diagrams. 

In this paper we study minimization diagrams for another 
class of convex functions, different from the convex distance 
functions. The functions we study take the form f(x,y) = 
g(x) + h(y), where g and h are triply-differentiable and 
g'(x)g"'(x) < g"{x) 2 and h'{y)ti"{y) < h"{y) 2 for all 
x and y. For example, squared Euclidean distance has this 
form with g(x) — x 2 and h{y) = y 2 . More generally, 
f(x,y) — \x\ c + \y\ c (for c > 1) satisfies these requirements Q 
and its minimization diagrams are the Voronoi diagrams fol- 
ic distance, known to have linear complexity (TJ. We show 
in Section [III] that the minimization diagrams of arbitrary 
functions in this class have analogous properties: 

• Any two level sets S r (p) = {q \ f(q — p) — r} are 
simple closed curves that intersect in at most two points; 
if they intersect at two points, they cross properly at these 
points. That is, these sets, which are defined analogously 

'The divergence of the double derivative of \x\ c at the origin when e < 2, 
and its non-differentiability there when 2 < c < 3, do not present any serious 
difficulties to our theory. 



to Euclidean circles, form a family of pseudocircles | 13 1, 

KH- 

Any bisector B(p, q) = {s \ f(s — p) = f(s — q)} forms 
either an axis-parallel line or a simple curve that is both 
a; -monotone and y-monotone; f(s — p) and f(s — q) vary 
unimodally along the bisector B(p 1 q). 
Any two bisectors B(p,q) and B(p,q') intersect in at 
most one point; if they intersect, they do so at a proper 
crossing. That is, if we fix p and let q vary, the bisectors 
B(p, q) form a weak pseudoline arrangement j3j, |6j, 
|7j, and the correspondence between q and B(p,q) can 
be viewed as a form of duality for these pseudolines. 
Each cell of the minimization diagram is simply con- 
nected, as it is a single cell in a pseudoline arrangement. 
Thus, the minimization diagram has linear complexity 
and, like Voronoi diagrams for convex distance functions 
and abstract Voronoi diagrams^] can be constructed in 
randomized expected time O(nlogn) (Theorem pj). 



minimization diagrams for the convex function 

(1 + e*) 2 



f(x, y) = m 



ln(l + cosy). 



In Section IV we provide examples showing that our restric- 
tions are necessary: minimization diagrams for more general 
convex functions do not in general have these properties. 

Our motivation for studying this type of minimization 
diagram comes from smoothed distance. Given a fixed point 
o in the plane, the smoothed distance or biotope transform 
metric d (p,q) = 2d(p,q)/(d(p,o) + d(q,o) + d(p,q)) [?], 
is a geometric analogue of the Jaccard similarity measure 
for clustering binary data. A maximal set of points with a 
given minimum smoothed distance will be distributed around 
o with exponentially decreasing density [2|, but as we show 
in Section |VI| the point spacing may be improved by using 
smoothed distance in Lloyd's algorithm [?], ]T7) , a continuous 
variant of X-means clustering that repeatedly moves each site 
to the centroid of its Voronoi cell. 

Smoothed distance is a monotonic function of the dilation 
(d(p,o) + d(o,q))/d(p,q), a measure of the quality of a star 
network as an approximation to the Euclidean distances among 
a set of points. For the maximum dilation pair of points (p, q) 
from a set P for a fixed center o, q is among the O(l) nearest 
neighbors to p either in Euclidean distance or in the sequence 
of points sorted by distance from o, and this result forms 
the basis of an efficient algorithm for finding the center o 
that minimizes the maximum dilation [8j . However, this can 
be simplified using smoothed distance: the maximum-dilation 
pair must be neighbors in the smoothed distance Voronoi 
diagram (Proposition [T|. 

Neither smoothed distance nor dilation are translates of 
convex functions. However, in Section [TT] we use complex 
logarithms to transform the plane and we perform a suitable 
monotonic transformation of smoothed distance, showing that 
Voronoi diagrams for smoothed distance are equivalent to 



2 We do not show that these diagrams satisfy all requirements of an abstract 
Voronoi diagram, as we do not bound the number of crossings between 
bisectors for unrelated pairs of points. 



As we require, this function is the sum of univariate functions 
g(x) and h(y), with g' g'" < {g") 2 . It is not true that h'ti" < 
{h") 2 for all y, but it is true for —tt/2 < y < tt/2, so smoothed 
distance Voronoi diagrams are well behaved whenever each 
Voronoi cell spans an angle of at most tt/2 on each side of 
its site with respect to o (Theorem [9j. 

II. Dilation, Smoothed Distance, and Logarithmic 
Transformation 

Two important measures of similarity between sets A and 
B are the Hamming distance dn(A,B) = \A A B\ (where 
A denotes the symmetric difference of sets) and the Jaccard 
distance [12], which weights the Hamming distance by the 
size of the union of the sets: 

T(A R\ = lAABl = d H {A,B) 

1 ' ' \AVB\ d H (A,®)+d H (B,(!>)~d H (A,B)- 

A monotone transformation of this distance lies in the interval 

[0,1]: 

d.,(A,B) = — =- 



2d H (A,B) 



J(A,B) 



2 d H (A,<b) + d H (B,<b)+d H (A,B) 



The same formula defining the (modified) Jaccard distance in 
terms of the Hamming distance can be used to derive a new 
metric from any given metric space (X, d) and fixed point 
o 6 X. Define the o- smoothed distance or biotope transform 
metric by the formula 

, / n = g) 

o[P ' q) d(p,o) + d(q,o) + d(p,qy 

This is a metric on X \ {o}, as can be shown using the 
tight span H), fTT) , the minimal Loo-like metric completion 
of a metric space. Any four points {o,p, q,r} C (X, d) may 
be embedded isometrically into a metric space formed by an 
axis-aligned rectangle of the L\ plane (possibly a degener- 
ate rectangle), together with four line segments (possibly of 
length zero) connecting the corners of this rectangle to the 
four points. When the four line segments all have length 
zero, the four points are in cyclic order (o, p, r, q) on 
the rectangle's corners, and the rectangle has aspect ratio 
a, then the triangle inequality for d is satisfied exactly: 
d (p,q) = 1 = ^pi + ^pi = d Q {p,r) + d (q,r). Increasing 
the length of the line segments connecting the four points to 
the rectangle or placing the points in a different cyclic order 
only strengthens the triangle inequality. 

For the points within a d-ball of radius e d(o, c) around a 
point c, the factor d(p, o)+d(q, o) + d(p, q) in the definition of 
d ranges from (1 — e)2d(o, c) to (l + 2e)2d(o, c), varying by a 
factor of approximately l+3e within this range as e approaches 
zero. Thus, closely spaced sets of points in the smoothed 
distance have distances that are approximately similar to their 
unsmoothed distances. As a metric space on I\ {o}, the 
smoothed distance d a is topologically equivalent to the metric 



induced by d on the same space: both distances have the same 
open sets and neighborhood structures. 

Smoothed distance d for the Euclidean plane is invariant 
with respect to rotations and scaling centered at o. These 
transformations may be expressed by representing point (x, y) 
as the complex number x + iy, with o = 0: if z is any 
complex number, then the product zp may be interpreted 
geometrically as rotating the point p by the angle ZzOl 
and scaling the rotated point by \z\. With this interpretation, 
d Q (zp, zq) = d a (p,q). 

As Figure [T] shows, 
smoothed-distance circles 
with smaller radius closely 
resemble Euclidean circles, 
while larger-radius smoothed- 
distance circles are not even 
convex. 

Smoothed distance is 
closely related to dilation, a 
measure of the quality of a 
graph as an approximation 
to a metric space [5|. The 
dilation of vertices p and q is 
the ratio of their distances in 
the graph and in the ambient 

space; the dilation of the graph is the largest dilation of any 
pair of vertices. For a star graph G with center o and all 
other vertices as leaves |8|, the dilation of pair (p,q) is 




Fig. 1. Concentric circles for 
smoothed distance in the Euclidean 
plane, with o = (0,0), p = (1,0), 
and radii 0.5, 0.6, 0.7, 0.8, and 0.9 



d{p,o) + d(o,q) 



d(p,q) d a (p,q) 
and the dilation of the whole star graph is 



1 



max 

p,gev(G) 



d(p, o) + d(o, q) 
d(p, q) 



max — r- — 1. 

P,geV(G) d {p, q) 



Because dilation is a monotonic transformation of smoothed 
distance, the point p that is nearest to a query point q in 
terms of smoothed distance is also the point that has the 
greatest dilation with q. The Voronoi diagram for smoothed 
distance partitions the plane into regions such that the region 
containing a query point q corresponds to the point p for 
which a path from q will have to make the greatest detour 
by passing through o instead of connecting directly. The 
adjacency relations between cells in this Voronoi diagram are 
also meaningful for dilation: 

Proposition 1: The points p and q defining the dilation of 
a star graph have adjacent cells in the Voronoi diagram for 
smoothed distance, using the star's leaves as Voronoi sites and 
its hub as the point o. 

Proof: Form the Voronoi diagram for V(G) \ {q} and 
then add q to form the Voronoi diagram of V(G). Let R be 
the Voronoi region of p prior to adding q. R must contain q, 
because otherwise some other point than p would define the 
greatest dilation with respect to q. After adding q, part of R 
becomes incorporated into the Voronoi region for q, while the 




Fig. 2. Concentric circles for the logarithmically transformed distance S with 
radii 0.5, 0.75, 0.9, 0.99, and 0.999. 



rest of R (in particular the point p itself) remains in the region 
of p. Thus, these two regions meet. ■ 
Thus, we would like to understand and compute Voronoi 
diagrams for smoothed distance. Smoothed distance is neither 
convex nor translation-invariant, but we can transform it into 
an equivalent distance that is, by interpreting polar coordinates 
in the plane for which we are computing smoothed distance 
as Cartesian coordinates in a transformed plane. Equivalently, 
using complex numbers (with o = 0), if complex number z has 
polar coordinates (r, 8), then Inz has Cartesian coordinates 
(lnr, 8). Define a distance 6 on the transformed complex 
number plane by the formula 

S(p,q)=d„^ t e"). 

This logarithmic transformation replaces the invariance of 
d a under complex multiplication by invariance of 6 under 
complex addition (that is, translation of the plane): 

S(p + z,q + z) =S(p,q). 

Figure [2] shows concentric circles for 5. 

We may calculate 6((x, y), (0, 0)) directly as a formula 
of x and y: it equals d (p,q), where o — (0,0), p = 
(e x cos y, e x siny), and q = (1,0). Thus, d(p,o) = e x , 
d(q, o) — 1, and 



d(p, q) = \J (e x cos y — l) 2 + (e x sin y) 2 ) 
= \J (e x ) 2 — 2e x cos y + 1. 

Therefore, the smoothed distance is 

2d(p, q) 



d (p,q) = 



d{p,o) +d(q,o) + d(p, q) 



2y/(e x ) 2 - 2e x cosy + 1 
e x + 1 + \/{e x ) 2 -2e x cosy + 1 ' 

It is convenient to monotonically transform this formula: 



d (p,q) 



1 



(e x ) 2 + 26^ + 1 
(e x ) 2 — 2e x cos y + 1 



and 



2 l 1 l/ \d {p,q) \ 



(cosy + l)e x 
(e x ) 2 + 2e x + 1' 




-6 -4 -2 2 4 6 -3 -2 -1 1 2 3 



Fig. 3. Graphs of the functions g(x) and h(y) used to represent the 
transformed value of the smoothed distance. 

Therefore, if we define 

, ^ ) - b {\{ 1 - 1/ {dr)- 1 )'))- 

it follows that 

(1 + e x f 

f(x,y) =ln — ln(cosy + l). 

Thus, we have represented smoothed distance as a monotonic 
transform of translates of f(x, y) = g(x) + h(y), where 

(1 + e x ) 2 

9{x) = In and h(y) = - ln(cosy + 1). 

e 

Graphs of these two functions are depicted in Figure [3] 

III. Minimization Diagrams of Convex Functions 

The transformation in Section [II] from a Voronoi diagram 
of smoothed distance to a minimization diagram of translates 
of f(x, y) — g(x) + h(y) motivates the more general study of 
minimization diagrams of this type of function. A univariate 
function / : R R is convex if the set of points 

{(x,y)\y>f(x)} 

on or above f(x) is a convex subset of the plane. In higher 
dimensions, a multivariate function / : M. d \— > R is convex if 
the set 

{(x ,Xx,X 2 , ...Xd) | f(x ,Xx,X 2 , ■ ..Xd-l) <Xd} 

is a convex subset of R rf+1 . A convex univariate function 
is strictly convex if there is no interval within which it is 
linear. For a doubly-differentiable univariate function g(x), 
convexity is equivalent to the inequality g" > and strict 
convexity is equivalent to the inequality g" > 0. A curve in 
the (x, y)-plane is x-monotone if every line parallel to the y 
axis intersects it in at most one point; intuitively, these curves 
are the graphs of functions from x to y. Symmetrically, a curve 
is y-monotone if every line parallel to the x axis intersects it 
in at most one point. These concepts should be distinguished 
from that of a function being monotonic ally increasing or 
strictly monotonically increasing: if x\ > Xq, and 4>(x) is 
monotonically increasing, then 4>(x\) > 4>(xq); if <j>(x) is 
strictly monotonically increasing, then 4>{x\) > <f>(xo). In a 
minimization diagram of translates of a function f(x,y), a 
bisector B(p,q) of two distinct points p and q is the locus 
of points with equal values of f p and f q : B(p, q) = {r | 
f p (r) = f q (r)}. The minimization diagram of the points {p, q} 



is formed by partitioning the plane into two cells along the 
bisector. 

Proposition 2: Let f(x,y) — g(x) + h(y), where g and 
h are strictly convex. Then any bisector B(p, q) is either an 
axis-parallel line or an x-monotone and y-monotone curve. 

Proof: If p and q have equal x-coordinates, then f p (r) = 
f q (r) whenever h(r y —p y ) = h(r y —q y )', since this condition 
is independent of x, the bisector must be the union of one or 
more lines parallel to the y axis. Otherwise, any line x = C 
parallel to the y axis intersects the bisector B(p,q) in the 
points r such that h(r y — p y ) — h(r y — q y ) — g(C — q x ) — 
g(C — p x ). By strict convexity, the function H(y) = h(r y — 
Py) — h(r y — q y ) is strictly monotonically increasing, so there 
is at most one such point and the bisector is x-monotone. 
Symmetrically, if p and q have the same y-coordinate, the 
bisector must be the union of one or more lines parallel to 
the x axis, and otherwise it is y-monotone. The only curves 
that are simultaneously x-monotone or a union of y-parallel 
lines, and y-monotone or a union of x-parallel lines, are the 
ones listed in the proposition: a single axis-parallel line, or a 
doubly-monotone curve. ■ 

In particular, each bisector must be a pseudoline (the image 
of a line under a homeomorphism of the plan^ji, because 
it is either itself a line or a monotonic curve that partitions 
the plane into two cells. In the next sequence of lemmas, we 
show that the level sets of the translates of / are pseudocodes 
(simple closed curves that cross each other at most twice) by 
comparing their curvatures at tangent points of the same slope. 
For convenience, when the argument of the functions g(x) and 
h(y) is specified in the context, the notations g, h, g', h' , etc., 
refer to the values g(x), h(y), £y(x), ^h(y), etc. 

Lemma 3: Let /(x,y) = g(x) + h(y) where g and h are 
convex and differentiable. Let p = (x p , y p ) be any point other 
than the global minimum of /, let f(p) = r, and let S r be the 
level set {q | f(q) = r}. Then the slope of the tangent to S r 
at p is —h'/g'. 

Proof: Since p is not the minimum, S r is a simple 
closed curve. Therefore there is a unique tangent line which 
is characterized by the fact that, at points q on the line at a 
small distance e from p, the difference between f(p) and f(q) 
vanishes to first order. A point q at a distance proportional to e 
on the line of slope —h'/g' can be found by adding {—h'e, g'e) 
to p; at this point, f(q) is f (p) + (- h'e) (g'+ 0(e)) + (g'e) (h'+ 
0(e)) = f(p) + 0(e 2 ), so the difference vanishes to first order 
as desired. ■ 

Next, suppose we move e units along the curve S r from p; 
to within first order, the point we move to is found by adding 
(e ~ h =, e , 9 =) to p. When we do so, g' changes by a 

y/h' 2 +g' 2 y/g' 2 +h' 2 

factor of 1 — e — % = + 0(e 2 ), and h' changes by a factor 

g'\/g' 2 +h' 2 

of 1 + e — 9 h + 0(e 2 ), so the slope h'/g' changes by 

h'y<g' 2 + h' 2 

a factor of 1 + e 9 h ^Lj=L ^ 9 + 0(e 2 ). We may view the 

y/g' 2 +h' 2 

3 This definition is from |2l) ; see for a comparison with other common 
definitions of pseudolines. 



term - — / + a J g_ a pp ear j n g j n this expression as a local 

^Jg' 2 +h' 2 

measure of the curvature of S r ; it is not rotation-invariant 
but can be used to compare the curvatures of two curves at 
points of equal tangent slope. Larger values of this term mean 
a smaller radius of curvature and lower values mean a greater 
radius. To show that the radius of curvature at a given slope 
increases as r increases, we will examine the behavior of this 
function as we increase r by a small quantity e while moving 
p along a curve that keeps the slope of the tangent to S r fixed. 

Lemma 4: Let f(x, y) = g(x) + h(y) where g and h are 
convex and twice-differentiable. At any point p for which 
g',h' 7^ 0, define the curve T p through p to consist of points 
with the same value of —h'/g' as at p. Then the tangent line 
to T p at p has slope jr/^r. 

Proof: Move p by a distance proportional to some small 
e along this line, by adding the vector (e^r,e^) to p. As p 
moves, g' and h! both change by factors of l + e g g/ ^, + 0(e 2 ). 
Both factors are equal to within first order, so the slope —h'/g' 
does not change to first order and we have identified the correct 
tangent line to T p . ■ 

Lemma 5: Let f(x, y) — g(x) + h(y) where g and h are 
convex and triply-differentiable, and where g"'g' < (g") 2 and 
h"'h' < (h") 2 for all x and y. Then, along any curve T p , 
the radius of curvature of the curves S r at the points where 
these curves are crossed by T p is a monotonically increasing 
function of f(p). 

Proof: We assume without loss of generality that g, h, 
g', and b! are all positive at p, for otherwise we may achieve 
these assumptions by translating the plane, reflecting it across 
one or both of the coordinate axes, or adding a constant to /, 
without changing the truth of the lemma. As in the previous 
lemma, we move at a distance proportional to e along T p , 
to within first order, by adding to p. And as in 

the previous lemma, this motion causes g' and h' to both 



change by a factor of 1 + e 



g"h" 

g'h' 



0(e ); this same factor 



also describes the change in \J g' 2 + h' 2 . Thus, to find the 
direction of change to the value 9 - /^jj^L ll that we are 

using to compare curvatures for a given tangent slope, it only 
remains to evaluate the change to g" and h". The double 
derivative g", and therefore the term h'g"/g' appearing in 
the numerator of our local curvature function, changes by a 
factor of 1 + + 0{e 2 ); if g"'g' < [g") 2 , this is smaller 



than the factor of 1 + e 



g"h" 

g'h' 



0(e 2 ) describing the change to 
the denominator of our local curvature function. The double 
derivative h", and therefore the term g'h" jh! appearing in the 
numerator of our local curvature function, changes by a factor 
of l + e<|V +0(e 2 ); again, if h"'h' < [h") 2 ), this is smaller 



than the factor of 1 + e 9 g ,^, + 0(e 2 ) describing the change to 
the denominator of our local curvature function. Thus, if the 



assumptions of the lemma are met, 



g'h" /h'+h'g" / g' 



decreases 



vV 2 +?>' 5 

and the radius of curvature increases as f(p) increases along 



T p . 



Due to the convexity of /, and the strict convexity of g 
and h, the level sets S r (p) = {q | f p (q) = r} are themselves 
convex; in particular, they are simple closed curves in the plane 
and, as the following proposition shows, they form a family 
of pseudocodes in the plane. 

Proposition 6: Let f(x, y) — g(x) + h(y) be a convex func- 
tion such that g and h are triply-differentiable, g"'g' < (g") 2 , 
and h"'hf < (h") 2 , and let r\ and r 2 be numbers and p\ and 
P2 be points such that (ri,pi) ^ {'T2,P'i)- Then the level sets 
S ri (pi) intersect in at most two points; if they intersect at two 
points, they cross properly at these points. 

Proof: We suppose that a pair that forms more than two 
points of intersection or that has two non-crossing intersections 
exists, and proceed to derive a contradiction. If there are 
exactly two points of intersection, one of which is not a 
crossing, then (because both are simple closed curves) both 
must be points of tangency; we may increase the radius of 
the inner level set slightly and form two level sets that cross 
four times, so we may assume without loss of generality that 
there are three or more points of intersection. If r\ =Ti, then 
again we may change one of the radii slightly while preserving 
the property of having more than two points of intersection, 
for this change of radii cannot remove any crossings and can 
be chosen to turn at least half of the points of tangency into 
pairs of crossings. And if some of the points of intersection are 
tangencies, we may increase or decrease r\ by a small amount 
and replace at least half of the tangencies by crossings without 
changing the property that there are more than two intersection 
points. Thus, we may assume that the two level sets intersect 
more than two times at proper crossings. 

Now, let T be the set of real numbers x > 1 such that 
S ri (pi) and S r2 (xp2 + (1 — x)pi) have more than two proper 
crossings; that is, we consider translating p 2 directly away 
from pi for as far as we can while preserving the overly large 
number of crossings between the two curves. Because both 
level sets are bounded, / is itself bounded; let t = sup T and 
let ps = tp2 + (1 — t)pi. In order for the translated level sets 
to have more than two crossings for x < t but only to have 
two crossings for x > t, the sets S ri (j>i) and S r2 (ps) must 
be tangent. However, this point of tangency cannot be one at 
which the two level sets cross, because our assumptions imply 
that both curves have curvature that varies continuously along 
the curves. It cannot be a tangency in which the curve with the 
smaller value of r, lies inside the curve with the larger value 
of Ti, because then for x > t the tangency would become 
two crossings and t would not be equal to sup T. It cannot be 
a tangency in which the two curves meet externally, because 
then by convexity for x < t the curves would have only two 
crossing points near the tangency and again t would not be 
equal to supT. And it cannot be a tangency in which the 
curve with the smaller value of lies outside the curve with 
the larger value of r-j, because that would violate Lemma [5] 
However, these exhaust the possible ways the two curves can 
be tangent. This contradiction completes the proof. ■ 

The next proposition states that (with the same assumptions 
as Lemma [5] and Proposition [6]) the bisectors B(p,q) and 



B(p,s) act like pseudolines: they meet in at most one point, 
and if they meet they cross properly. 

Proposition 7: Let f(x, y) = g(x) + h(y) be a convex func- 
tion such that g and h are triply-differentiable, g"'g' < (g") 2 , 
and h"'h' < (h") 2 . Then any two bisectors B(p,q) and 
B(p, s) defined from / have at most one point of intersection. 
If they intersect, they cross properly. 

Proof: We show that B(p, q) and B(p, s) meet in at most 
one point by showing that any two points t and u in the plane 
are intersected by at most one bisector B(p,-). If there is a 
bisector B(p, q) that contains both of these points, then p and 
q are equidistant from t and from s; that is, p and q both 
belong to the level sets St = {s | fit — s) = f(t — p)} and 
Su = { s I f( u ~ s ) — f( u ~ P)}- These level sets are rotated 
by 7r from the level sets of Proposition [6] but that rotation does 
not affect the conclusion of the proposition: they have at most 
two points of intersection. One of these intersection points is 
p and the other is q; there can be no third intersection to form 
another bisector with p through t and u. If B(p, q) and B(p, s) 
met in two points, it would violate the uniqueness of bisectors 
through pairs of points, so such a double intersection cannot 
happen. 

To complete the proof of the proposition, we must show 
that if two bisectors intersect, they meet in a proper crossing. 
But if two bisectors B(p,q) and B(p,s) met at a point of 
tangency without crossing, then a small translation of s either 
towards or away from p would transform this tangency into a 
pair of crossings, violating the first part of the proposition. ■ 

In other words, if we fix p and let q vary, the family 
of bisectors B(p, q) forms a weak pseudoline arrangement. 
However, the proposition applies only to pairs of bisectors that 
share one of their defining points. These results are not quite 
enough to show that minimization diagrams for / are abstract 
Voronoi diagrams in the sense of Klein [14], because abstract 
Voronoi diagrams require all bisectors to have a constant 
number of intersection points. However, the same general 
results as for abstract Voronoi diagrams follow in this case. 

Theorem 8: Let f(x, y) — g(x)+h(y) be a convex function 
such that g and h are triply-differentiable, g"'g' < (g") 2 , and 
h"'h' < (h") 2 . Then any minimization diagram for /, defined 
by a finite set P of n point sites, subdivides the plane into 
n simply-connected regions, one per site. The diagram can 
be constructed in time 0(n log n) using a primitive that finds 
minimization diagrams for three sites. 

Proof: We may assume without loss of generality that 
the minima of g and h occur at x = y = 0, for otherwise 
we may add an appropriate constant to x and y without 
changing the combinatorial description of the minimization 
diagram. Thus, the cell for any point p contains p itself. In any 
weak arrangement of pseudolines, any nonempty intersection 
of halfspaces defined by the pseudolines forms a single cell of 
the arrangement (e.g. see Theorem 11.4.11 of (7J). It follows 
by Proposition [7] that the cell of p in the minimization diagram 
is a single cell in a weak arrangement of pseudolines and is 
therefore simply connected. Thus, the minimization diagram 
for / and P consists of n simply connected cells. 



We construct the diagram by a standard randomized incre- 
mental algorithm in which we add points to P one at a time 
according to a uniformly random permutation, and maintain a 
history DAG describing the cells of the diagram in past states 
of the construction. We also maintain the sequence of bisectors 
surrounding each cell of the diagram, in a balanced binary tree 
data structure. To add a point p, we use the history DAG to 
locate the cell of the diagram that contains p. We then form 
a list L of cells known to overlap with the cell of p; initially 
this list includes only the cell containing p. To build the cell 
for p, we repeatedly remove the cell for a site q from L, and 
split this cell along the bisector B(p,q). The points where 
this bisector crosses the boundary of the cell for q may be 
found by a binary search of the boundary, in which each step 
consists of finding the vertex of the minimization for p, q, and 
a third site determining one of the boundary segments, and 
comparing the x and y coordinates of this vertex to those of 
the other vertices on the segment. After the part to be removed 
from the cell of q is determined, the boundary segments in the 
binary search tree for q are removed and the associated cells 
are added to L if necessary. 

The time to locate p is O(logri) in expectation by a 
standard analysis of history DAGs. Each new feature of the 
minimization diagram takes O(logn) time to construct, using 
the binary search trees, and there are in expectation 0(1) new 
features for the ith added site. Thus, the total expected time 
for the construction is 0(n log n). ■ 

IV. Two Bad Examples 

The assumptions that we make on the form of f(x,y) = 
g(x) + h(y) as a sum of two univariate convex functions, 
and on the triple derivatives of these functions, may seem 
technical and unnecessary. However, in this section we provide 
examples showing that without the assumption on the form 
of / beyond convexity, its minimization diagrams may have 
quadratic complexity, and without assumption on the triple 
derivatives of g and h, the level sets for / may not be 
pseudocircles. 

Figure [4] shows an 
example in which g(x) 
and h(y) grow so quickly 
that g"'g' > {g") 2 
and ti"ti > (ti 1 ) 2 
for sufficiently large 
x. Specifically, 
//"// .:l(i.r' • 
while [g") 2 




(16a; 4 + 16a; 2 + 4)e 2x ; it 
follows that g"'g' > {g") 2 
for x > 1/V2. One level 
set has been translated so 
that it crosses the outer 
level set four times, so 
these level sets are not 

pseudocircles. More generally, whenever g"'g' > (g") 2 
over some interval of values of x, the level sets of 



Fig. 4. The level sets for fix, y) = 
2 2 

e x + e v (shown for function values 
2.5, 5, 10, 20, 40, 80, and 160) do not 
form pseudocircles. 



f(x,y) = g(x) + g(y) do not form pseudocircles: the radius 
of curvature at one of the tangents with slope —1 shrinks 
rather than growing for increasing values of /. 

Figures [5] and [6] provide a sketch of a construction for a 
convex function that has minimization diagrams with quadrat- 
ically growing complexity. The function grows most slowly on 
a horizontal line through the origin, and more quickly along 
any other horizontal line, so that for any sites p and q that 
are not on a horizontal line, the Voronoi cell for p dominates 
that for q for points far enough away from both sites on a 
horizontal line through p. In particular the vertical line of 
cells to the right of the figure generates a sequence of cells 
with horizontal boundaries that extends, despite interruptions, 
through the whole figure. A more widely spaced sequence of 
sites at the bottom of the figure interrupts these cells, splitting 
them into linearly many pieces. 




Fig. 5. Level sets for a convex function the minimization diagrams of which 
have quadratic complexity. 




Fig. 6. Sketch of the structure for a minimization diagram of a convex 
function with quadratic complexity. 

V. Voronoi Diagrams for Smoothed Distance 

So, now that we've transformed the smoothed distance 
Voronoi diagram into the form of a minimization diagram for 
a function f(x, y) = g(x) + h(y), and now that we know con- 
ditions on g and h that ensure that the minimization diagram 
to have linear complexity and be constructable efficiently, do 
the specific g and h arising in this application meet these 
conditions? The answer is: it depends. 

First, consider the function g(x) = \n(e x + 2 + e~ x ), with 
derivatives g'{x) = {e x -l)/(e x + 1), g"(x) = 2e x /(l + e x ) 2 , 
andg"'(x) = —2e x (e x — l)/(e x + l) 3 . Forpositive x, g'" and 
therefore g"'g' are negative while (g") 2 is positive, so g"'g' < 
(g") 2 - Because the function is symmetric (g(— x) = g{x)) the 
same holds for negative x. And at zero, g"'g' — while (g") 2 
is nonzero. Therefore, g meets the conditions of Theorem [8] 

Next, consider the function h(y) = — ln(l + cosy). Its 
derivatives are h'(y) — tan(y/2) = siny/(l + cosy), 
h"{y) = 1/(1 + cosy), and ti"(y) = sin y/(l + cosy) 2 . Then 
ti"h' = sin 2 y/(l+cosy) 3 , while {h") 2 = l/(l+cosy) 2 . The 
ratio of these two values, sin 2 y/(l + cosy) 2 = tan 2 (y/2), 
is less than one when |y| < tt/2 and greater than one 
when \y\ > ir/2. Thus, h satisfies the requirement that 
h"'h' < (h") 2 only for y in the interval (-tt/2, tt/2). This 
implies that we may only apply Theorem [8] to smoothed 



distance when each Voronoi cell consists of points spanning 
at most a right angle with o and the cell's site. 

Theorem 9: Let P be a point set such that, in the Voronoi 
diagram for o-smoothed distance, every point q within the 
Voronoi cell for a site p forms an angle Zpoq that is at most 
tt/2. Then each cell in the Voronoi diagram is connected, and 
the diagram may be constructed in randomized expected time 
0(n\ogn). 

Proof: To prove this, in outline, we replace the point set 
P by its logarithmically transformed image {p | e p € P}; each 
point in P corresponds to infinitely many transformed points, 
all with the same x coordinate and with y coordinates that 
differ by integer multiples of 2tt. Under this transformation, 
each point in the transformed plane is associated with a 
Voronoi region that the exponential function maps back to the 
correct Voronoi region for the associated input point in the 
original plane. Then, instead of using the functions g(x) and 
h(y) in the transformed plane, as defined above, we replace the 
function h(y) by a modified function that has the same values 
within the interval (—tt/2, tt/2) but that obeys the inequality 
h"'h' < (h") 2 for larger values as well. This replacement 
allows us to apply Theorem [8] to the transformed input, and 
we show that, with the assumptions stated in the theorem, it 
produces the same cell decomposition as the one we wish 
to compute. We use an efficient randomized incremental 
algorithm to compute the smoothed Voronoi diagram, similar 
to the algorithm used in Theorem [8] 

Smoothed distance may be replaced by the translates of 
a convex function by logarithmically mapping the sites and 
monotonically transforming the distance values. Thus, Voronoi 
diagrams for smoothed distance are closely related to min- 
imization diagrams for this convex function. The specific 
relation is this: if we view the points in the plane as complex 
numbers, with o = 0, and map the finite set P of sites 
to the infinite vertically-periodic set L = {q \ e q £ P}, 
then the exponential function forms a covering map from 
the minimization diagram of / with respect to L to the 
Voronoi diagram for smoothed distance of P. The cells in the 
minimization diagram are mapped many-to-one to cells in the 
Voronoi diagram, and edges and features in the minimization 
diagram are mapped to edges and features in the Voronoi 
diagram, etc. Although it is problematic to perform geometric 
algorithms on infinite point sets such as L, we may use our 
analysis to show that cells in the minimization diagram for 
L are simply connected, while performing the randomized 
incremental algorithm described in the proof of Theorem [8] 
directly on the finite point set P. 

There is a difficulty with this approach, however: our proof 
that the cells are simply connected does not apply, because 
the convex function into which we have transformed smoothed 
distance does not meets the requirement of Theorem [8] that the 
function h satisfies the inequality h"'h' < (h") 2 , at least for 
large values of y. And although the assumptions of Theorem [9] 
imply that only small values of y need be considered in the 
final Voronoi diagram, larger values may be needed at early 
stages of our randomized incremental construction. 



Therefore, rather than computing the Voronoi diagram for 
smoothed distance d Q (p,q) itself, we compute the Voronoi 
diagram for a (non-metric) distance function D(p, q) = 
g(ln(d(p,o)/d(q,o))) + h(Zpoq). Here g(x) = \n(e x + 2 + 
e~ x ), matching the monotonically-transformed formula for 
smoothed distance. However, we set h(y) to be a function that 
equals — ln(l + cosy) for \y\ < ir/2 but that is a quadratic 
polynomial for values of y outside this range; the two quadratic 
polynomials (for positive and negative y) are chosen to have 
a value, first derivative, and second derivative that matches 
— ln(l + cosy) at — 7r/2 and 7r/2. However, being quadratic 
polynomials, they have third derivative equal to zero, and 
therefore satisfy the requirement that h"'h' < (h") 2 for the 
range of values of y in which they define the value of h. 
The fact that the overall h function has a discontinuous third 
derivative at —ir/2 and at ir/2 does not cause any problems 
for our analysis. 

With this modified distance function, Theorem [8] can be 
used to show that the cells of the minimization diagram 
for g(x) + h(y) with respect to the points of L are simply 
connected, and therefore that their images, the cells of the 
Voronoi diagram for D with respect to the points of P, are 
also connected. We may then apply a randomized incremental 
algorithm of the type described in Theorem [8] to construct 
the Voronoi diagram for D with respect to the points of P. 
At intermediate stages of this construction, the cells of the 
diagram may have self-adjacencies (corresponding to bound- 
aries between two different images in I of a point in P) but 
these need not be treated any differently than any other of the 
bisectors in the diagram, and must vanish before the algorithm 
finishes. 

It remains to show that each cell in the diagram for the 
modified distance function is equal to the corresponding cell 
in the diagram for the true smoothed distance d a . Observe that, 
for any site p, the cell for p for distance D is contained within 
the boundaries of the cell for p for distance d : by assumption, 
these boundaries are all within the region of the plane in which 
D(p,q) and (the corresponding monotonic transformation of) 
d (p,q) coincide, so on those boundary points d (p,q) is 
accurately represented by D(p, q) while the distances from 
q to other points as measured by D may be underestimates 
of the true value as measured by d Q . Thus, the cells for the 
Voronoi diagram for D are subsets of the corresponding cells 
for the Voronoi diagram of d D . However, since the cells for 
the Voronoi diagram for D nevertheless cover the plane, both 
sets of cells coincide and the computed Voronoi diagram is 
correct. ■ 

VI. Lloyd's Algorithm 

Evenly spaced points in Euclidean and related metric spaces 
have applications ranging from coding theory (T7J and color 
quantization flO) to dithering (spatial halftoning) and stip- 
pling for image rendering fl8) , (22). However, random points 
typically have uneven spacing, and metric e-nets (maximal 
point sets such that no two points are closer than e to each 
other), while more uniform, may still have varying density. 



A common method for improving the spacing of Euclidean 
point sets is Lloyd's algorithm [17], a variant of the K- 
means clustering algorithm that repeatedly computes a Voronoi 
diagram and replaces each point by the centroid of its Voronoi 
cell. Its output is a centroidal Voronoi diagram [?], a well- 
spaced collection of points that form the centroids of their 
Voronoi cells. 

Clarkson [i2| suggested using metric e-nets for o-smoothed 
distance to generate well-spaced point sets with a distribution 
centered at o that decreases exponentially with distance from 
o. For this application, one must restrict the points to an 
annulus centered on o; otherwise, an e-net would not have a 
finite number of points. As an alternative means of generating 
exponentially-distributed and well-spaced points in this annu- 
lus, we experimented with a variant of Lloyd's algorithm that 
uses smoothed distance in place of Euclidean distance in its 
calculations. Specifically, rather than computing a Euclidean 
Voronoi diagram of the given points, we computed a Voronoi 
diagram for the smoothed distance. And rather than moving 
each point to the centroid of its cell C (the point p minimizing 
I q ec d(P> ( ?) 2< ^C)' we move each point to the point minimizing 
/ eC d (p,q) 2 dC. Finally, in order to make the measure dC 
of area used in the definition of the area integral transform 
scale-invariantly to match the symmetries of the smoothed 
distance, we chose a measure that is uniform not in the 
Euclidean plane in which the smoothed distance is defined, 
but rather the uniform measure in the transformed plane that 
has the Euclidean polar coordinates (log r, 8) as its Cartesian 
coordinates. 

For simplicity, our implementation performs its calculations 
in the Euclidean plane, rasterized as a bitmap image. We 
compute the Voronoi diagram by finding the nearest site 
to each pixel of the rasterized annulus, and approximate 
f qec l) 2<1C h y S geC d o{Pi <l) 2 /d(q, o) 2 , where the sum 
is over the pixels in the Voronoi region C and the l/d(q,o) 2 
term weights each pixel by its measure in the transformed 
plane. The results of one run of our implementation are 
depicted in Figure [7] We found that the iteration quickly 
(within two iterations) smoothed out any gross variation in the 
spacing of the given points, and then more slowly converged 
to a more ideal shape for each Voronoi cell. It was necessary 
to choose an initial set of points that was exponentially 
distributed around o; we placed each point by choosing L 
uniformly within the interval [In r, In R] (where r and R are 
the inner and outer radius of the annulus) and 8 uniformly in 
[0, 2n], and then selecting the point (e L cos8, e L sin#) in the 
Euclidean plane. We found that if instead we selected points 
with uniform Euclidean measure in the annulus, too many 
points were placed far from o and too few were placed close 
to o; Lloyd's algorithm was slow in correcting this imbalance. 

VII. Conclusions 

We have identified a general condition on convex func- 
tions that causes their minimization diagrams to have linear 
complexity, applied this condition to Voronoi diagrams of 
smoothed distance and to finding the minimum dilation pair 





Fig. 7. Lloyd's algorithm on an annulus with an 18:1 ratio of inner to 
outer radius, for 128 exponentially-distributed random points. Top: initial 
configuration and one iteration. Center: after two and four iterations. Bottom: 
after 4, 8, and 16 iterations. The large black dots are the sites at each iteration; 
the smaller white spots represent the smoothed-distance Voronoi centroids to 
which these sites will be moved at the next iteration. 



of leaves in a star network, and experimented with using 
a smoothing algorithm based on these Voronoi diagrams to 
generate evenly-spaced points exponentially distributed around 
a given center point. 

Several directions for further research remain open: 

• If we translate the convex function f(x, y) in three 
dimensions rather than two by adding independent con- 
stants to the values for each point site, when does the 
resulting planar minimization diagram still have linear 
complexity? For instance, additively weighting f(x, y) = 
x 2 + y 2 in this way results in a power diagram. For 
any convex f(x,y) = g(x) + h(y), the bisectors of 
an additively weighted minimization diagram are still 
monotone curves, but we can no longer guarantee that 
they form pseudoline arrangements. Nevertheless it might 
be possible that they form minimization diagrams with 
connected cells. 

• Is it possible to characterize the functions that can 
be monotonically transformed into the form f(x, y) = 
g(x) + h(y) with g and h both convex? The functions that 



have this form already are exactly the convex functions 
for which every axis-parallel rectangle has equal sums on 
its two pairs of opposite corners. However, if a function 
does not already have this form it may not be clear how 
to transform it into this form, as we did for smoothed 
distance. 

. Does our condition g"'g' < (g") 2 and h"'h' < (h") 2 
characterize the convex functions g and h such that the 
translated level sets of f(x,y) — g(x) + h(y) form 
pseudocodes? It does when g = h: If g"'g' > (g") 2 , 
then level sets of f(x, y) = g(x) + g(y) have a curvature 
at slope 1 that grows tighter for larger circles. However 
the situation is less clear when g and h may differ. 

• In particular, does the convex function f(x,y) = g(x) + 
h(y) coming from smoothed distance and dilation have 
level sets that are pseudocircles when y > 7r/2? 

• If not, is there some other natural distance function 
d on the nonzero complex numbers, satisfying scale 
invariance d(zp, zq) — d(j>, q), that has simply connected 
Voronoi regions for all sets of two or more points in 
general position? Note that the most obvious choice, 
the Euclidean distance between bap and bag, does not 
work: in general there will be a single Voronoi region 
containing the origin, which will not be simply connected. 
The replacement function used in the proof of Theorem [9] 
does not seem very natural. 

• If we define a maximization diagram in which the sites 
are point pairs (p, q) and the function to be maximized 
is the dilation d(p, q) /(d(p, o) + d(q, o)) of p and q with 
respect to a query point o, our previous results [ 8 1 imply 
that this diagram has 0(n) cells. Are these cells simply 
connected? 

• Can we characterize the convex functions whose level sets 
are pseudocircles? For instance, as well as the functions 
g(x) + h(y) studied here, this is also true of convex 
distance functions. 

• If a convex function has level sets that are pseudocircles, 
do its minimization diagrams automatically have simply 
connected cells, or is an additional condition required for 
this to be true? 

• To what extent can this theory be generalized to mini- 
mization diagrams in higher dimensions? 
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